%
% @brief 射频系统的 冲激响应特性 补偿
% @author 41204@qq.com / 2025-10
% @note 本类中的函数均以 Static 方式定义，以方便使用。
%

classdef hcomp
    methods(Static, Access = public)
        %% 补偿射频接收机的幅度和相位失真
        % 输入：
        %   A: 输入IQ信号 (复数)
        %   B: 输出IQ信号 (复数)
        % 输出：
        %   H_estimated: 估计的系统频率响应
        %   A_compensated: 补偿后的信号
        function [H_estimated, A_compensated] = compensate_rf_response(A, B, varargin)        
            p = inputParser;
            addParameter(p, 'smoothing_factor', 0.01, @(x) x>=0 && x<=1);
            addParameter(p, 'compensation_limit_db', 20, @isnumeric);
            parse(p, varargin{:});
            
            N = length(A);
            
            % 1. 估计系统频率响应
            A_fft = fft(A, N);
            B_fft = fft(B, N);
            
            % 计算原始频率响应
            H_raw = B_fft ./ A_fft;
            
            % 处理除零和噪声问题
            small_values = abs(A_fft) < max(abs(A_fft)) * 1e-6;
            H_raw(small_values) = 1;  % 在这些频率点设为无失真
            
            % 2. 平滑频率响应（重要！）
            H_smoothed = rm.rf_hfreq_compensate.smooth_frequency_response(H_raw, p.Results.smoothing_factor);
            H_estimated = H_smoothed;
            
            % 3. 设计补偿滤波器（逆系统）
            H_comp = 1 ./ H_estimated;
            
            % 限制最大补偿幅度（避免噪声放大）
            max_comp_mag = 10^(p.Results.compensation_limit_db/20);
            H_comp_mag = abs(H_comp);
            H_comp(H_comp_mag > max_comp_mag) = max_comp_mag .* exp(1j * angle(H_comp(H_comp_mag > max_comp_mag)));
            
            % 4. 应用补偿
            B_fft_comp = B_fft .* H_comp;
            A_compensated = ifft(B_fft_comp, 'symmetric');
        end
    end

    methods (Static, Access = private)
        %% 使用移动平均或低通滤波平滑频率响应
        function H_smoothed = smooth_frequency_response(H_raw, smoothing_factor)
            if smoothing_factor > 0
                window_size = max(3, round(length(H_raw) * smoothing_factor));
                if mod(window_size, 2) == 0
                    window_size = window_size + 1;
                end
                
                % 幅度平滑
                mag_smoothed = movmean(abs(H_raw), window_size);
                
                % 相位平滑（需要小心处理相位卷绕）
                phase_raw = angle(H_raw);
                phase_unwrapped = unwrap(phase_raw);
                phase_smoothed = movmean(phase_unwrapped, window_size);
                
                H_smoothed = mag_smoothed .* exp(1j * phase_smoothed);
            else
                H_smoothed = H_raw;
            end
        end
    end
end
